% Plot the forcasted seismicity across the Aftershock Zone of May 2008
% earthquake in SW Iceland

ni= 1 ; %: length(vec_m)
forecast_data = N_robust_p50_xy(:,:,ni); vec_m(ni) 

% Prob of earthquake for the entire aftershock zone
Plb = P_mean ; 
 
min(min(forecast_data))
max(max(forecast_data))
caxis_max = 0.05; 

%% 
Dxx = [0.0,-0.1,0.1,-0.1,0.1,-0.1];   
Dyy = [-0.15,0.10,-0.1,0.10,-0.1,0.10]; 
dxx = [0.02,0.01,0.01,0.01,0.01,0.01];
dyy = [0.00,0.01,0.01,0.01,0.01,0.01];
X_upperText = lonMin-2*0.02;% 13.22;
Y_upperText = latMax+3*0.01;% 42.50-0.05/2;

%% Plot Forecasts N
gcf=figure;
set(gcf,'Position',[200   80   900   700])

pcolor(Xcgrid,Ycgrid,forecast_data);
labels = ['Forecasted Seismicity for Earthquakes with M', num2str(vec_m(ni)),'^{+}'];
hcb = colorbar; 
caxis([0.0 caxis_max])
cmap = colormap(hot(256)); % jet(5)
cmap = flipud(cmap); 
colormap(cmap)
ylabel(hcb,labels,'FontSize',21,'Rotation',-270,'FontWeight','normal'); hcb.Label.Position(1) = 4;
shading interp
hold on

% Plot Regions
Regions = shaperead('SIce_shp/S.IceTowns_prj.shp'); % geographic data structure array in projected map coordinates
vec_Regions   = 1:42; 
for j=1:length(vec_Regions)
    jj = vec_Regions(j);
    plot(Regions(jj).X,Regions(jj).Y,'b','linewidth',1.3)
    hold on
end

deltaT_figure = 1;
deltaT_figure_scale = '[day]';

for k= 2:length(vec_m)
    if Plb(k) >= 0.10
        text(lonMin+2*deltaGrid_X,latMax-.02-0.03*(k-2),['Pr(M\geq',num2str(vec_m(k)),') = ',num2str(Plb(k),'%2.2f')],'FontName','Times New Roman','FontSize',18,'FontWeight','normal')
    elseif (Plb(k) >= 0.010 && Plb(k) < 0.10)
        text(lonMin+2*deltaGrid_X,latMax-.02-0.03*(k-2),['Pr(M\geq',num2str(vec_m(k)),') = ',num2str(Plb(k),'%3.2f')],'FontName','Times New Roman','FontSize',18,'FontWeight','normal')
    else
        text(lonMin+2*deltaGrid_X,latMax-.02-0.03*(k-2),['Pr(M\geq',num2str(vec_m(k)),') = ',num2str(Plb(k),'%10.1e')],'FontName','Times New Roman','FontSize',18,'FontWeight','normal')
    end
end

h1=rectangle('Position',[lonMin,latMin+.01,lonMax-lonMin-.05,latMax-latMin-.01],'FaceColor','none','EdgeColor',[.5 .5 .5],'LineWidth',3.5);

text(X_upperText,Y_upperText,{['{\it T}_{start}= ',num2str(0,'%4.2f'),'UTC ',...
  datestr(addtodate(datenum( '05/29/2008','mm/dd/yy'),1,'day'),'dd-mm-yy')...
  ' FI= ',num2str(deltaT_figure,3),' ',deltaT_figure_scale]},'FontName','Times New Roman','FontSize',17,'FontWeight','normal')

xlabel('Longitude','fontsize',24,'FontWeight','normal')
ylabel('Latitude','fontsize',24,'FontWeight','normal')
xlim([lonMin-0.1 lonMax+.1])
ylim([latMin-0.02 latMax+0.07])

set(gca,'fontsize',20)
set(gca,'XTick',str2double(num2str(min(xlim),'%4.1f')):0.2:str2double(num2str(max(xlim),'%4.1f')),...
    'YTick',str2double(num2str(min(ylim),'%4.1f')):0.1:str2double(num2str(max(ylim),'%4.1f')),'GridLineStyle',':')

saveas(gcf,[pwd,'/','forecast_Nm2-N50Pmean.tif']) 



